### Replication files for
### "The Local Costs of Moving: How Residential Moves Weaken Local Engagement"
### Hans Lueders
### Comparative Political Studies
### September 2025




# Set-up ------------------------------------------------------------------

### packages
library(ggplot2)
library(foreign)
require(lmtest)
library(stargazer)
library(xtable)
library(maptools)
library(gridExtra)
library(ggmap)
library(rgdal)
library(sp)
library(gridExtra)
library(grid)
library(lattice)
library(Matching)
library(ggpubr)
library(geosphere)
library(spatstat)
library(raster)
library(lfe)
library(AER)
library(dplyr)
library(corrplot)
library(ggrepel)
library(rbounds)
library(mgcv)
library(mediation)
library(sensemakr)
library(cjoint)
library(emmeans)
library(cregg)
library(mapproj)
library(readxl)
library(sf)
library(Hmisc)
library(estimatr)
library(lme4)



### empty environment
rm(list=ls())


### set working directory
path <- "/Users/hanslueders/Dropbox/Publications/14 CPS_The Local Costs of Moving/Replication files/" # define path here
setwd(path)






# Load and prepare datasets -----------------------------------------------

### Load data
load("data/TheLocalCostsOfMoving_ReplicationData.RData")






# Figure 1: Changes in engagement before and after a move -----------------

### plot
p <- ggplot() + 
  geom_hline(yintercept = 0, lty=2, size=0.8) +
  geom_vline(xintercept = -0.5, lty=2, size=0.8) +
  geom_vline(xintercept = 0.5, lty=2, size=0.8) +
  geom_point(data=soep_fig1, 
             aes(x=years, y=b, col = years), size=2.6,
             position = position_dodge(width=0.75)) +
  geom_errorbar(data=soep_fig1, 
                aes(x=years, ymin=(b-1.96*se), ymax=(b+1.96*se), col = years), 
                position = position_dodge(width=0.75), width=0, linewidth=1.5) +
  facet_wrap(~outcomeFACTOR, scales="free_y") +
  scale_x_continuous(breaks=seq(-10,10,2)) +
  scale_color_gradient2(low = "#c6393c", high = "#4747d1", mid="gray70") +    
  xlab("Years before/after move") + 
  ylab("Change in the difference\nMovers vs. stayers") +  
  theme_bw() +
  theme(plot.title = element_text(size=18, face="bold", hjust=.5),
        axis.title = element_text(size=17, face="bold"),
        axis.text = element_text(size=15),
        strip.text = element_text(size=17, face="bold"),
        legend.position = "none")
p
ggsave("Figure1_ChangeEngagement_main.png", p, width=14,height=7)




# Figure 2: Changes in engagement around moves of varying distances -------

### plot
p <- ggplot() + 
  annotate("rect", xmin=2.5,xmax=3.5, ymin=-Inf,ymax=Inf, alpha=.5, fill="seashell3") +
  geom_hline(yintercept = 0, lty=2, size=0.8) +
  geom_point(data=soep_fig2, 
             aes(x=as.numeric(yearsFACTOR), y=b, col=modelFACTOR, shape=modelFACTOR), size=3,
             position = position_dodge(width=0.5)) +
  geom_errorbar(data=soep_fig2, 
                aes(x=as.numeric(yearsFACTOR), ymin=(b-1.96*se), ymax=(b+1.96*se), col=modelFACTOR), 
                position = position_dodge(width=0.5), width=0, size=1.5) +
  scale_x_continuous(breaks=1:5,
                     labels = c("baseline","1 to 5\nyears\nbefore", "year\nof\nmove",
                                "1 to 5\nyears\nafter","6+\nyears\nafter")) +
  scale_color_manual(values = c("gray70", "gray60", "gray40", "gray20"), name="") +
  scale_shape_manual(values = c(15,19,17,8), name="") +
  facet_wrap(~outcomeFACTOR, scales="free_y") +
  ylab("Change in the difference\nbetween\nmovers and stayers") +  
  theme_bw() +
  theme(plot.title = element_text(size=18, face="bold", hjust=.5),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=17, face="bold"),
        axis.text = element_text(size=15),
        strip.text = element_text(size=17, face="bold"),
        legend.title = element_blank(),
        legend.text = element_text(size=15),
        legend.position = "bottom")
p
ggsave("Figure2_ChangeEngagement_distance.png", p, width=14,height=7)





# Figure 3: Most moves occur over short distances -------------------------

### Plots
# Distances moved
grav2015_mig$title <- "Distance"
p1 <- ggplot() + 
  geom_density(data = grav2015_mig, aes(x=distance/1000),
               fill = "gray70", col=NA, alpha=.85) +
  geom_vline(xintercept = median(grav2015_mig$distance,na.rm=T)/1000, lty=2, size=1) +  
  xlab("Distance moved (kms)\n(all moves)") + 
  ylab("Density") +
  scale_x_continuous(limits=c(0,800), breaks=seq(0,800,100)) +
  facet_wrap(~title) +  
  theme_bw() +
  theme(strip.text = element_text(size=15, face="bold"),
        axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.position = "none") 


# Share moving to neighboring county
dist.agg <- aggregate(grav2015_mig$bothContiguous, by=list(grav2015_mig$orig_kreis), FUN=mean, na.rm=T)
colnames(dist.agg) <- c("AGS", "shareContiguous")
dist.agg$title <- "Neighboring counties"
p2 <- ggplot() + 
  geom_density(data=dist.agg, aes(x=shareContiguous*100), fill="gray70", col=NA, alpha=0.85) + 
  geom_vline(xintercept = median(dist.agg$shareContiguous,na.rm=T)*100, lty=2, size=1) +    
  scale_x_continuous(limits=c(0,100), breaks=seq(0,100,20)) +
  facet_wrap(~title) +
  xlab("Share of all moves (%)\n(by county)") + 
  ylab("Density") + 
  theme_bw() +
  theme(strip.text = element_text(size=15, face="bold"),
        axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.position = "none") 

# Share moving within the same state
grav2015_mig$orig_land <- NA
grav2015_mig$orig_land[grav2015_mig$orig_kreis <= 9999] <- substr(as.character(grav2015_mig$orig_kreis[grav2015_mig$orig_kreis <= 9999]), 1,1)
grav2015_mig$orig_land[grav2015_mig$orig_kreis >= 10000] <- substr(as.character(grav2015_mig$orig_kreis[grav2015_mig$orig_kreis >= 10000]), 1,2)

grav2015_mig$dest_land <- NA
grav2015_mig$dest_land[grav2015_mig$dest_kreis <= 9999] <- substr(as.character(grav2015_mig$dest_kreis[grav2015_mig$dest_kreis <= 9999]), 1,1)
grav2015_mig$dest_land[grav2015_mig$dest_kreis >= 10000] <- substr(as.character(grav2015_mig$dest_kreis[grav2015_mig$dest_kreis >= 10000]), 1,2)

grav2015_mig$sameland <- ifelse(grav2015_mig$orig_land == grav2015_mig$dest_land, 1, 0)

dist.agg <- aggregate(grav2015_mig$sameland, by=list(grav2015_mig$orig_kreis), FUN=mean, na.rm=T)

dist.agg$title <- "Same state"
p3 <- ggplot() + 
  geom_density(data=dist.agg, aes(x=x*100), fill="gray70", col=NA, alpha=0.85) + 
  geom_vline(xintercept = median(dist.agg$x,na.rm=T)*100, lty=2, size=1) +    
  scale_x_continuous(limits=c(0,100), breaks=seq(0,100,20)) +
  facet_wrap(~title) +
  xlab("Share of all moves (%)\n(by county)") + 
  ylab("Density") + 
  theme_bw() +
  theme(strip.text = element_text(size=15, face="bold"),
        axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        legend.position = "none") 


# combining all three plots
p <- ggarrange(p1,p2,p3,nrow=1)
p
ggsave("Figure3_ShortDistances.png", p, width=12, height=3.5)





# Figure 4: Movers tend to move between politically similar enviro --------

##### Figure 4a: national engagement #####
### Create data to be plotted
# calculate difference in turnout
grav2015_mig$BT_turnout_diff <- abs(grav2015_mig$orig_BT_turnout - grav2015_mig$dest_BT_turnout)
grav2015_pop$BT_turnout_diff <- abs(grav2015_pop$orig_BT_turnout - grav2015_pop$dest_BT_turnout)


# Define BT election result differences
grav2015_mig$BT_second_diff <- 0.5*(abs(grav2015_mig$orig_BT_second_union_share - grav2015_mig$dest_BT_second_union_share) + 
                                         abs(grav2015_mig$orig_BT_second_spd_share - grav2015_mig$dest_BT_second_spd_share) + 
                                         abs(grav2015_mig$orig_BT_second_fdp_share - grav2015_mig$dest_BT_second_fdp_share) + 
                                         abs(grav2015_mig$orig_BT_second_greens_share - grav2015_mig$dest_BT_second_greens_share) +
                                         abs(grav2015_mig$orig_BT_second_linke_share - grav2015_mig$dest_BT_second_linke_share) + 
                                         abs(grav2015_mig$orig_BT_second_afd_share - grav2015_mig$dest_BT_second_afd_share))

grav2015_pop$BT_second_diff <- 0.5*(abs(grav2015_pop$orig_BT_second_union_share - grav2015_pop$dest_BT_second_union_share) + 
                                         abs(grav2015_pop$orig_BT_second_spd_share - grav2015_pop$dest_BT_second_spd_share) + 
                                         abs(grav2015_pop$orig_BT_second_fdp_share - grav2015_pop$dest_BT_second_fdp_share) + 
                                         abs(grav2015_pop$orig_BT_second_greens_share - grav2015_pop$dest_BT_second_greens_share) +
                                         abs(grav2015_pop$orig_BT_second_linke_share - grav2015_pop$dest_BT_second_linke_share) + 
                                         abs(grav2015_pop$orig_BT_second_afd_share - grav2015_pop$dest_BT_second_afd_share))


# create datasets
dmig1 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, BT_turnout_diff))
dpop1 <- subset(grav2015_pop, select = c(orig_kreis, dest_kreis, BT_turnout_diff))
dmig1$sample <- 1
dmig1$median <- median(dmig1$BT_turnout_diff, na.rm=T)
dpop1$sample <- 2
dpop1$median <- median(dpop1$BT_turnout_diff, na.rm=T)
d1 <- rbind(dmig1, dpop1)
d1$sample <- factor(d1$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d1$outcome <- 1
colnames(d1)[3] <- "diff"

dmig2 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, BT_second_diff))
dpop2 <- subset(grav2015_pop, select = c(orig_kreis, dest_kreis, BT_second_diff))
dmig2$sample <- 1
dmig2$median <- median(dmig2$BT_second_diff, na.rm=T)
dpop2$sample <- 2
dpop2$median <- median(dpop2$BT_second_diff, na.rm=T)
d2 <- rbind(dmig2, dpop2)
d2$sample <- factor(d2$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d2$outcome <- 2
colnames(d2)[3] <- "diff"

d <- rbind(d1, d2)
d$outcome <- factor(d$outcome, levels = 1:2,
                    labels = c("Federal election turnout",
                               "Federal election result"))



### Plot
p <- ggplot() + 
  geom_density(data=d, aes(x=diff, fill=sample), col=NA, alpha=.7) + 
  geom_vline(data=d, aes(xintercept = median, lty=sample), size=1) +
  facet_wrap(~outcome, scales="free") +
  scale_x_continuous(breaks = seq(0,0.5,0.05)) +
  scale_fill_manual(values = c("steelblue4", "gray70"), name="") +
  scale_linetype_manual(values = c(2,3), name="") +
  xlab("Difference between origin and destination") + 
  ylab("Density") +
  theme_bw() +
  theme(axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        strip.text = element_text(size=15, face="bold"),
        legend.position = "bottom",
        legend.text = element_text(size=12))   
p
ggsave("Figure4a_differences_national.png", p, width=10,height=4)




##### Figure 4b: local engagement #####

### Create data to be plotted
# Define difference in LT election turnout
grav2015_mig$LT_turnout_diff <- abs(grav2015_mig$orig_LT_turnout_lag - grav2015_mig$dest_LT_turnout_lag)
grav2015_pop$LT_turnout_diff <- abs(grav2015_pop$orig_LT_turnout_lag - grav2015_pop$dest_LT_turnout_lag)

# Define difference in LT election results
grav2015_mig$LT_diff <- 0.5*(abs(grav2015_mig$orig_LT_votes_union_share_lag - grav2015_mig$dest_LT_votes_union_share_lag) + 
                                  abs(grav2015_mig$orig_LT_votes_spd_share_lag - grav2015_mig$dest_LT_votes_spd_share_lag) + 
                                  abs(grav2015_mig$orig_LT_votes_fdp_share_lag - grav2015_mig$dest_LT_votes_fdp_share_lag) + 
                                  abs(grav2015_mig$orig_LT_votes_greens_share_lag - grav2015_mig$dest_LT_votes_greens_share_lag) +
                                  abs(grav2015_mig$orig_LT_votes_linke_share_lag - grav2015_mig$dest_LT_votes_linke_share_lag) + 
                                  abs(grav2015_mig$orig_LT_votes_afd_share_lag - grav2015_mig$dest_LT_votes_afd_share_lag))

grav2015_pop$LT_diff <- 0.5*(abs(grav2015_pop$orig_LT_votes_union_share_lag - grav2015_pop$dest_LT_votes_union_share_lag) + 
                                  abs(grav2015_pop$orig_LT_votes_spd_share_lag - grav2015_pop$dest_LT_votes_spd_share_lag) + 
                                  abs(grav2015_pop$orig_LT_votes_fdp_share_lag - grav2015_pop$dest_LT_votes_fdp_share_lag) + 
                                  abs(grav2015_pop$orig_LT_votes_greens_share_lag - grav2015_pop$dest_LT_votes_greens_share_lag) +
                                  abs(grav2015_pop$orig_LT_votes_linke_share_lag - grav2015_pop$dest_LT_votes_linke_share_lag) + 
                                  abs(grav2015_pop$orig_LT_votes_afd_share_lag - grav2015_pop$dest_LT_votes_afd_share_lag))

# Create difference in CSO density
grav2015_mig$CSO_diff <- abs(grav2015_mig$orig_CSO_stock_count_pc - grav2015_mig$dest_CSO_stock_count_pc)
grav2015_pop$CSO_diff <- abs(grav2015_pop$orig_CSO_stock_count_pc - grav2015_pop$dest_CSO_stock_count_pc)



# create datasets
dmig1 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, LT_turnout_diff))
dpop1 <- subset(grav2015_pop, select = c(orig_kreis, dest_kreis, LT_turnout_diff))
dmig1$sample <- 1
dmig1$median <- median(dmig1$LT_turnout_diff, na.rm=T)
dpop1$sample <- 2
dpop1$median <- median(dpop1$LT_turnout_diff, na.rm=T)
d1 <- rbind(dmig1, dpop1)
d1$sample <- factor(d1$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d1$outcome <- 1
colnames(d1)[3] <- "diff"

dmig2 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, LT_diff))
dpop2 <- subset(grav2015_pop, select = c(orig_kreis, dest_kreis, LT_diff))
dmig2$sample <- 1
dmig2$median <- median(dmig2$LT_diff, na.rm=T)
dpop2$sample <- 2
dpop2$median <- median(dpop2$LT_diff, na.rm=T)
d2 <- rbind(dmig2, dpop2)
d2$sample <- factor(d2$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d2$outcome <- 2
colnames(d2)[3] <- "diff"

dmig3 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, CSO_diff))
dpop3 <- subset(grav2015_pop, select = c(orig_kreis, dest_kreis, CSO_diff))
dmig3$sample <- 1
dmig3$median <- median(dmig3$CSO_diff, na.rm=T)
dpop3$sample <- 2
dpop3$median <- median(dpop3$CSO_diff, na.rm=T)
d3 <- rbind(dmig3, dpop3)
d3$sample <- factor(d3$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d3$outcome <- 3
colnames(d3)[3] <- "diff"

d <- rbind(d1,d2,d3)
d$outcome <- factor(d$outcome, levels = 1:3,
                    labels = c("State election turnout",
                               "State election result",
                               "Civil society organizations"))


### Plot
p <- ggplot() + 
  geom_density(data=d, aes(x=diff, fill=sample), col=NA, alpha=.7) + 
  geom_vline(data=d, aes(xintercept = median, lty=sample), size=1) +
  facet_wrap(~outcome, scales="free") +
  scale_fill_manual(values = c("steelblue4", "gray70"), name="") +
  scale_linetype_manual(values = c(2,3), name="") +
  xlab("Difference between origin and destination") + 
  ylab("Density") +
  theme_bw() +
  theme(axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        strip.text = element_text(size=15, face="bold"),
        legend.position = "bottom",
        legend.text = element_text(size=12))   
p
ggsave("Figure4b_differences_local.png", p, width=12,height=4)





# Table 1: Political differences are associated with smaller bilat --------

### Regressions
grav$diff <- grav$diff_BT_turnout
m1 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)

grav$diff <- grav$diff_BT_result
m2 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)

grav$diff <- grav$diff_LT_turnout
m3 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)

grav$diff <- grav$diff_LT_result
m4 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)

grav$diff <- grav$diff_CSO
m5 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)


### Stargazer output
stargazer(m1,m2,m3,m4,m5,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          column.labels = c("Federal turnout", "Federal results", "State turnout", "State results", "CSOs"),
          dep.var.labels.include = T,
          dep.var.labels  = c("Emigration (log)"),
          covariate.labels=c("Difference in engagement",
                             "Distance (log kms)",
                             "Diff. manufacturing",
                             "Diff. services",
                             "Diff. GDP p.c.",
                             "Diff. UER",
                             "Diff. patents",                             
                             "GDP (log) at origin",
                             "GDP (log) at destination",
                             "Population (log) at origin",
                             "Population (log) at destination"),           
          digits = 3,
          add.lines = list(
            c("Origin-FE?", rep("YES", 5)),
            c("Destination-FE?", rep("YES", 5)),
            c("Year-FEs?", "YES", rep("YES", 5))),
          notes = c("Standard errors clustered by dyad"))





# Table 2: Smaller geographic scale ---------------------------------------

### Regressions
# Difference in BT turnout
gravmun$diff <- gravmun$diff_BT_turnout
m1 <- felm(count_log ~ diff +
             distance_km_log +
             orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
           data=gravmun)
m2 <- felm(count_log ~ diff +
             distance_km_log +
             orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
           data=gravmun,
           subset=gravmun$samekreis == 1)

# Difference in BT result
gravmun$diff <- gravmun$diff_BT_result
m3 <- felm(count_log ~ diff +
             distance_km_log +
             orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
           data=gravmun)
m4 <- felm(count_log ~ diff +
             distance_km_log +
             orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
           data=gravmun,
           subset=gravmun$samekreis == 1)

# Difference in LT turnout
gravmun$diff <- gravmun$diff_LT_turnout
m5 <- felm(count_log ~ diff +
             distance_km_log +
             orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
           data=gravmun)
m6 <- felm(count_log ~ diff +
             distance_km_log +
             orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
           data=gravmun,
           subset=gravmun$samekreis == 1)

# Difference in LT result
gravmun$diff <- gravmun$diff_LT_result
m7 <- felm(count_log ~ diff +
             distance_km_log +
             orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
           data=gravmun)
m8 <- felm(count_log ~ diff +
             distance_km_log +
             orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
           data=gravmun,
           subset=gravmun$samekreis == 1)

# Difference in CSOs
gravmun$diff <- gravmun$diff_CSO
m9 <- felm(count_log ~ diff +
             distance_km_log +
             orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
           data=gravmun)
m10 <- felm(count_log ~ diff +
              distance_km_log +
              orig_pop_total_log + dest_pop_total_log | orig_gemeinde + dest_gemeinde + year | 0 | dyad_FE2, 
            data=gravmun,
            subset=gravmun$samekreis == 1)


### stargazer output
stargazer(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          column.labels = c("Federal turnout", "Federal results", "State turnout", "State results", "CSOs"),
          column.separate = rep(2,5),
          dep.var.labels.include = T,
          dep.var.labels  = c("Emigration (log)"),
          covariate.labels=c("Difference in engagement",
                             "Distance (log kms)",
                             "Pop. at origin (log)",
                             "Pop. at dest. (log)"),           
          digits = 3,
          add.lines = list(
            c("Sample", rep(c("all", "local"), 5)),
            c("Origin-FE?", rep("YES", 10)),
            c("Destination-FE?", rep("YES", 10)),
            c("Year-FEs?", rep("YES", 10))),
          notes = c("Standard errors clustered by dyad"))






# Table A1: Descriptive statistics ----------------------------------------

### See STATA do-file





# Table A2: Relationship between domestic migration and survey dro --------

### See STATA do-file





# Figure A1: Distribution of measures of engagement -----------------------

### get variables to plot
# federal turnout
d1 <- subset(dcounty, year == 2021, select = c(BT_turnout))
colnames(d1) <- c("value")
d1$value <- d1$value * 100
d1$variable <- 1

# state turnout
d2 <- subset(dcounty, year == 2021, select = c(LT_turnout_lag))
colnames(d2) <- c("value")
d2$value <- d2$value * 100
d2$variable <- 2

# CSO stock
d3 <- subset(dcounty, year == 2021, select = c(CSO_stock_count_pc))
colnames(d3) <- c("value")
d3$variable <- 3


# combine
d.plot <- rbind(d1,d2,d3)
d.plot$variable <- factor(d.plot$variable, levels = 1:3,
                          labels = c("Federal election turnout", 
                                     "State election turnout",
                                     "CSOs per 1,000 capita"))


### plot
p <- ggplot() + 
  geom_density(data = d.plot[d.plot$variable != "Local CSOs per 1,000 capita",], 
               aes(x=value),
               fill = "gray70", col=NA, alpha=.85) +
  xlab("") + 
  ylab("Density") +
  facet_wrap(~variable, scales="free") +  
  theme_bw() +
  theme(strip.text = element_text(size=13, face="bold"),
        axis.title = element_text(size=12, face="bold"),
        axis.text = element_text(size=11),
        legend.position = "none") 
p
ggsave("FigureA1_distribution_engagement.png", p, width=9, height=4)





# Figure A2: Maps of measures of engagement -------------------------------

### load map shapefile
# shapefile
shp <- st_read("data/VG5000_KRS.shp") %>%  
  st_transform(4326)

# shapefile for state borders
shp.state <- st_read("data/VG5000_LAN.shp") %>%
  st_transform(4326)
shp.state <- subset(shp.state, GF == 9)


### add net migration data to shapefile
sub <- subset(dcounty, year == 2021, select=c(AGS, BT_turnout, 
                                        LT_turnout_lag, 
                                        CSO_stock_count_pc))
shp2 <- merge(shp, sub, by="AGS", all.x=T)


### prepare coordinates for labels 
centroids.df <- data.frame(st_geometry(st_centroid(shp2)), shp2$AGS)
colnames(centroids.df) <- c("geometry", "AGS")
centroids.df$name <- NA
centroids.df$name[centroids.df$AGS == "02000"] <- "Hamburg"
centroids.df$name[centroids.df$AGS == "03241"] <- "Hanover"
centroids.df$name[centroids.df$AGS == "05314"] <- "Cologne"
centroids.df$name[centroids.df$AGS == "06412"] <- "Frankfurt"
centroids.df$name[centroids.df$AGS == "08111"] <- "Stuttgart"
centroids.df$name[centroids.df$AGS == "09162"] <- "Munich"
centroids.df$name[centroids.df$AGS == "09564"] <- "Nuremberg"
centroids.df$name[centroids.df$AGS == "16051"] <- "Erfurt"
centroids.df$name[centroids.df$AGS == "14612"] <- "Dresden"
centroids.df$name[centroids.df$AGS == "14713"] <- "Leipzig"
centroids.df$name[centroids.df$AGS == "11000"] <- "Berlin"
centroids.df$name[centroids.df$AGS == "13003"] <- "Rostock"
centroids.df <- subset(centroids.df, !is.na(centroids.df$name))



### Maps
# Federal turnout
shp2$title <- "Federal election turnout"
p1 <- ggplot() + 
  geom_sf(data=shp2, aes(fill=BT_turnout*100), linewidth=.35) +
  geom_sf(data=shp.state, fill=NA, col="black", linewidth=.75) +
  scale_fill_gradient(low="white", high="black", na.value="gray90", name="Federal turnout") +
  geom_sf(data=centroids.df$geometry, fill="firebrick3", shape=24, size=2) +
  geom_sf_label(data=centroids.df$geometry, label=centroids.df$name, nudge_y = .25, size=5) +
  facet_wrap(~title) +  
  theme_bw() + 
  theme(axis.title=element_blank(), 
        axis.text=element_blank(), 
        axis.ticks=element_blank(),
        axis.line=element_blank(),
        strip.text = element_text(size=20, face="bold"),        
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        plot.title = element_text(size=20, face="bold", hjust=.5),
        panel.background = element_rect(fill = 'gray95', colour = NA),
        legend.position="bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=17),
        legend.key.width = unit(4, "line"))  

# State turnout
shp2$title <- "State election turnout"
p2 <- ggplot() + 
  geom_sf(data=shp2, aes(fill=LT_turnout_lag*100), linewidth=.35) +
  geom_sf(data=shp.state, fill=NA, col="black", linewidth=.75) +
  scale_fill_gradient(low="white", high="black", na.value="gray90", name="State turnout") +
  geom_sf(data=centroids.df$geometry, fill="firebrick3", shape=24, size=2) +
  geom_sf_label(data=centroids.df$geometry, label=centroids.df$name, nudge_y = .25, size=5) +
  facet_wrap(~title) +
  theme_bw() + 
  theme(axis.title=element_blank(), 
        axis.text=element_blank(), 
        axis.ticks=element_blank(),
        axis.line=element_blank(),
        strip.text = element_text(size=20, face="bold"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        plot.title = element_text(size=20, face="bold", hjust=.5),
        panel.background = element_rect(fill = 'gray95', colour = NA),
        legend.position="bottom",
        legend.title = element_blank(), 
        legend.text = element_text(size=17),
        legend.key.width = unit(4, "line"))  

# CSOs per capita
shp2$title <- "CSOs per 1,000 capita"
p3 <- ggplot() + 
  geom_sf(data=shp2, aes(fill=CSO_stock_count_pc), linewidth=.35) +
  geom_sf(data=shp.state, fill=NA, col="black", linewidth=.75) +
  scale_fill_gradient(low="white", high="black", na.value="gray90", name="CSOs per capita") +
  geom_sf(data=centroids.df$geometry, fill="firebrick3", shape=24, size=2) +
  geom_sf_label(data=centroids.df$geometry, label=centroids.df$name, nudge_y = .25, size=5) +
  facet_wrap(~title) +
  theme_bw() + 
  theme(axis.title=element_blank(), 
        axis.text=element_blank(), 
        axis.ticks=element_blank(),
        axis.line=element_blank(),
        strip.text = element_text(size=20, face="bold"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        plot.title = element_text(size=20, face="bold", hjust=.5),
        panel.background = element_rect(fill = 'gray95', colour = NA),
        legend.position="bottom",
        legend.title = element_blank(), 
        legend.text = element_text(size=17),
        legend.key.width = unit(4, "line"))  

### combine to one plot
p <- ggarrange(p1,p2,p3, nrow=1)
p
ggsave("FigureA2_maps_engagement.png", p, width=14, height=8)






# Table A3: Changes in engagement around a move, by distance of th --------

### See STATA do-file





# Table A4: Engagement after a move, by age  ------------------------------

### See STATA do-file






# Figure A3: Changes in engagement around a move, by education ------------

### plot
p <- ggplot() + 
  annotate("rect", xmin=2.5,xmax=3.5, ymin=-Inf,ymax=Inf, alpha=.5, fill="seashell3") +
  geom_hline(yintercept = 0, lty=2, size=0.8) +
  geom_point(data=soep_figA3, 
             aes(x=as.numeric(yearsFACTOR), y=b, col=educ_level, shape=educ_level), size=3,
             position = position_dodge(width=0.6)) +
  geom_errorbar(data=soep_figA3, 
                aes(x=as.numeric(yearsFACTOR), ymin=(b-1.96*se), ymax=(b+1.96*se), col=educ_level), 
                position = position_dodge(width=0.6), width=0, size=1.5) +
  scale_x_continuous(breaks=1:5,
                     labels = c("baseline","1 to 5\nyears\nbefore", "year\nof\nmove",
                                "1 to 5\nyears\nafter","6+\nyears\nafter")) +
  scale_color_manual(values = rep(c("gray60", "gray20"), each=2), name="") +
  scale_shape_manual(values = rep(c(15,19),2), name="") +
  facet_wrap(~outcomeFACTOR, scales="free_y") +
  ylab("Change in the difference\nbetween\nmovers and stayers") +  
  theme_bw() +
  theme(plot.title = element_text(size=18, face="bold", hjust=.5),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=17, face="bold"),
        axis.text = element_text(size=15),
        strip.text = element_text(size=17, face="bold"),
        legend.title = element_blank(),
        legend.text = element_text(size=15),
        legend.position = "bottom")
p
ggsave("FigureA3_ChangeEngagement_education.png", p, width=14,height=7)






# Table A5: Changes in engagement around a move, by level of educa --------

### See STATA do-file





# Figure A4: Movers tend to move between politically similar envir --------

##### Figure A4a: national engagement #####

### Create data to be plotted
# define difference in federal election turnout
grav2015_alt$BT_turnout_diff <- abs(grav2015_alt$orig_BT_turnout - grav2015_alt$dest_BT_turnout)

# define difference in federal election result
grav2015_alt$BT_second_diff <- 0.5*(abs(grav2015_alt$orig_BT_second_union_share - grav2015_alt$dest_BT_second_union_share) + 
                                      abs(grav2015_alt$orig_BT_second_spd_share - grav2015_alt$dest_BT_second_spd_share) + 
                                      abs(grav2015_alt$orig_BT_second_fdp_share - grav2015_alt$dest_BT_second_fdp_share) + 
                                      abs(grav2015_alt$orig_BT_second_greens_share - grav2015_alt$dest_BT_second_greens_share) +
                                      abs(grav2015_alt$orig_BT_second_linke_share - grav2015_alt$dest_BT_second_linke_share) + 
                                      abs(grav2015_alt$orig_BT_second_afd_share - grav2015_alt$dest_BT_second_afd_share))

# create datasets
dmig1 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, BT_turnout_diff))
dalt1 <- subset(grav2015_alt, select = c(orig_kreis, dest_kreis, BT_turnout_diff))
dmig1$sample <- 1
dmig1$median <- median(dmig1$BT_turnout_diff, na.rm=T)
dalt1$sample <- 2
dalt1$median <- median(dalt1$BT_turnout_diff, na.rm=T)
d1 <- rbind(dmig1, dalt1)
d1$sample <- factor(d1$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d1$outcome <- 1
colnames(d1)[3] <- "diff"

dmig2 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, BT_second_diff))
dalt2 <- subset(grav2015_alt, select = c(orig_kreis, dest_kreis, BT_second_diff))
dmig2$sample <- 1
dmig2$median <- median(dmig2$BT_second_diff, na.rm=T)
dalt2$sample <- 2
dalt2$median <- median(dalt2$BT_second_diff, na.rm=T)
d2 <- rbind(dmig2, dalt2)
d2$sample <- factor(d2$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d2$outcome <- 2
colnames(d2)[3] <- "diff"

d <- rbind(d1, d2)
d$outcome <- factor(d$outcome, levels = 1:2,
                    labels = c("Federal election turnout",
                               "Federal election result"))



### Plot
p <- ggplot() + 
  geom_density(data=d, aes(x=diff, fill=sample), col=NA, alpha=.7) + 
  geom_vline(data=d, aes(xintercept = median, lty=sample), size=1) +
  facet_wrap(~outcome, scales="free") +
  scale_x_continuous(breaks = seq(0,0.5,0.05)) +
  scale_fill_manual(values = c("steelblue4", "gray70"), name="") +
  scale_linetype_manual(values = c(2,3), name="") +
  xlab("Difference between origin and destination") + 
  ylab("Density") +
  theme_bw() +
  theme(axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        strip.text = element_text(size=15, face="bold"),
        legend.position = "bottom",
        legend.text = element_text(size=12))   
p
ggsave("FigureA4a_differences_national_alt.png", p, width=10,height=4)





##### Figure A4b: national engagement #####

### Create data to be plotted
# define difference in state election turnout
grav2015_alt$LT_turnout_diff <- abs(grav2015_alt$orig_LT_turnout_lag - grav2015_alt$dest_LT_turnout_lag)

# define difference in state election result
grav2015_alt$LT_diff <- 0.5*(abs(grav2015_alt$orig_LT_votes_union_share_lag - grav2015_alt$dest_LT_votes_union_share_lag) + 
                               abs(grav2015_alt$orig_LT_votes_spd_share_lag - grav2015_alt$dest_LT_votes_spd_share_lag) + 
                               abs(grav2015_alt$orig_LT_votes_fdp_share_lag - grav2015_alt$dest_LT_votes_fdp_share_lag) + 
                               abs(grav2015_alt$orig_LT_votes_greens_share_lag - grav2015_alt$dest_LT_votes_greens_share_lag) +
                               abs(grav2015_alt$orig_LT_votes_linke_share_lag - grav2015_alt$dest_LT_votes_linke_share_lag) + 
                               abs(grav2015_alt$orig_LT_votes_afd_share_lag - grav2015_alt$dest_LT_votes_afd_share_lag))

# define difference in CSO density
grav2015_alt$CSO_diff <- abs(grav2015_alt$orig_CSO_stock_count_pc - grav2015_alt$dest_CSO_stock_count_pc)




# create datasets
dmig1 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, LT_turnout_diff))
dalt1 <- subset(grav2015_alt, select = c(orig_kreis, dest_kreis, LT_turnout_diff))
dmig1$sample <- 1
dmig1$median <- median(dmig1$LT_turnout_diff, na.rm=T)
dalt1$sample <- 2
dalt1$median <- median(dalt1$LT_turnout_diff, na.rm=T)
d1 <- rbind(dmig1, dalt1)
d1$sample <- factor(d1$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d1$outcome <- 1
colnames(d1)[3] <- "diff"

dmig2 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, LT_diff))
dalt2 <- subset(grav2015_alt, select = c(orig_kreis, dest_kreis, LT_diff))
dmig2$sample <- 1
dmig2$median <- median(dmig2$LT_diff, na.rm=T)
dalt2$sample <- 2
dalt2$median <- median(dalt2$LT_diff, na.rm=T)
d2 <- rbind(dmig2, dalt2)
d2$sample <- factor(d2$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d2$outcome <- 2
colnames(d2)[3] <- "diff"

dmig3 <- subset(grav2015_mig, select = c(orig_kreis, dest_kreis, CSO_diff))
dalt3 <- subset(grav2015_alt, select = c(orig_kreis, dest_kreis, CSO_diff))
dmig3$sample <- 1
dmig3$median <- median(dmig3$CSO_diff, na.rm=T)
dalt3$sample <- 2
dalt3$median <- median(dalt3$CSO_diff, na.rm=T)
d3 <- rbind(dmig3, dalt3)
d3$sample <- factor(d3$sample, levels = 1:2,
                    labels = c("actual (dashed)", "benchmark (dotted)"))
d3$outcome <- 3
colnames(d3)[3] <- "diff"

d <- rbind(d1,d2,d3)
d$outcome <- factor(d$outcome, levels = 1:3,
                    labels = c("State election turnout",
                               "State election result",
                               "Civil society organizations"))


### Plot
p <- ggplot() + 
  geom_density(data=d, aes(x=diff, fill=sample), col=NA, alpha=.7) + 
  geom_vline(data=d, aes(xintercept = median, lty=sample), size=1) +
  facet_wrap(~outcome, scales="free") +
  scale_fill_manual(values = c("steelblue4", "gray70"), name="") +
  scale_linetype_manual(values = c(2,3), name="") +
  xlab("Difference between origin and destination") + 
  ylab("Density") +
  theme_bw() +
  theme(axis.title = element_text(size=14, face="bold"),
        axis.text = element_text(size=12),
        strip.text = element_text(size=15, face="bold"),
        legend.position = "bottom",
        legend.text = element_text(size=12))   
p
ggsave("FigureA4b_differences_local_alt.png", p, width=12,height=4)






# Table A6: gravity model -- cross-section --------------------------------

### Regressions
grav$diff <- grav$diff_BT_turnout
m1 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, subset=grav$year %in% c(2013))

grav$diff <- grav$diff_BT_result
m2 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, subset=grav$year %in% c(2013))

grav$diff <- grav$diff_LT_turnout
m3 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, subset=grav$year %in% c(2013))

grav$diff <- grav$diff_LT_result
m4 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, subset=grav$year %in% c(2013))

grav$diff <- grav$diff_CSO
m5 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, subset=grav$year %in% c(2013))


### Stargazer output
stargazer(m1,m2,m3,m4,m5,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          column.labels = c("Federal turnout", "Federal results", "State turnout", "State results", "CSOs"),
          dep.var.labels.include = T,
          dep.var.labels  = c("Emigration (log)"),
          covariate.labels=c("Difference in engagement",
                             "Distance (log kms)",
                             "Diff. manufacturing",
                             "Diff. services",
                             "Diff. GDP p.c.",
                             "Diff. UER",
                             "Diff. patents",                             
                             "GDP (log) at origin",
                             "GDP (log) at destination",
                             "Population (log) at origin",
                             "Population (log) at destination"),           
          digits = 3,
          add.lines = list(
            c("Origin-FE?", rep("YES", 5)),
            c("Destination-FE?", rep("YES", 5)),
            c("Year-FEs?", "YES", rep("YES", 5))),
          notes = c("Standard errors clustered by dyad"))





# Table A7: gravity model -- within- vs. cross-state moves ----------------

### Regressions
grav$diff <- grav$diff_BT_turnout
m1 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
           subset=grav$samestate == 1)
m2 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
           subset=grav$samestate == 0)

grav$diff <- grav$diff_BT_result
m3 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
           subset=grav$samestate == 1)
m4 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
           subset=grav$samestate == 0)

grav$diff <- grav$diff_LT_turnout
m5 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
           subset=grav$samestate == 1)
m6 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
           subset=grav$samestate == 0)

grav$diff <- grav$diff_LT_result
m7 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
           subset=grav$samestate == 1)
m8 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
           subset=grav$samestate == 0)

grav$diff <- grav$diff_CSO
m9 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
           subset=grav$samestate == 1)
m10 <- felm(count_emig2_log ~ diff + 
              distanceKM_log + 
              diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
              orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
              orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav, 
            subset=grav$samestate == 0)


### Stargazer output
stargazer(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          column.labels = c("Federal turnout", "Federal results", "State turnout", "State results", "CSOs"),
          column.separate = rep(2,5),
          dep.var.labels.include = T,
          dep.var.labels  = c("Emigration (log)"),
          covariate.labels=c("Difference in engagement",
                             "Distance (log kms)",
                             "Diff. manufacturing",
                             "Diff. services",
                             "Diff. GDP p.c.",
                             "Diff. UER",
                             "Diff. patents",                             
                             "GDP (log) at origin",
                             "GDP (log) at destination",
                             "Population (log) at origin",
                             "Population (log) at destination"),           
          digits = 3,
          add.lines = list(
            c("Sample", rep(c("same state", "diff. state"), 5)),
            c("Origin-FE?", rep("YES", 10)),
            c("Destination-FE?", rep("YES", 10)),
            c("Year-FEs?", "YES", rep("YES", 10))),
          notes = c("Standard errors clustered by dyad"))




# Table A8: gravity model - federal election turnout ----------------------

### Regressions
grav$diff <- grav$diff_BT_turnout
m1 <- felm(count_emig2_log ~ diff + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m2 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m3 <- felm(count_emig2_log ~ diff + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m4 <- felm(count_emig2_log ~ diff + 
             distanceKM_log +  
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m5 <- felm(count_emig2_log ~ diff + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m6 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)


### Output
stargazer(m1,m2,m3,m4,m5,m6,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          column.labels = c("Federal turnout"),
          column.separate = 6,
          dep.var.labels.include = T,
          dep.var.labels  = c("Emigration (log)"),
          covariate.labels=c("Difference",
                             "Distance (log kms)",
                             "Diff. manufacturing",
                             "Diff. services",
                             "Diff. GDP p.c.",
                             "Diff. UER",
                             "Diff. patents",
                             "Share manufacturing at origin",
                             "Share manufacturing at destination",
                             "Share services at origin",
                             "Share services at destination",
                             "GDP pc at origin",
                             "GDP pc at destination",
                             "UER at origin",
                             "UER at destination",    
                             "Patents pc at origin",
                             "Patents pc at destination",   
                             "GDP (log) at origin",
                             "GDP (log) at destination",
                             "Population (log) at origin",
                             "Population (log) at destination"),           
          digits = 3,
          add.lines = list(
            c("Origin-FE?", rep("YES", 6)),
            c("Destination-FE?", rep("YES", 6)),
            c("Year-FEs?", "YES", rep("YES", 6))),
          notes = c("Standard errors clustered by dyad"))




# Table A9: gravity model -- federal election result ----------------------

### Regressions
grav$diff <- grav$diff_BT_result
m1 <- felm(count_emig2_log ~ diff + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m2 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m3 <- felm(count_emig2_log ~ diff + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m4 <- felm(count_emig2_log ~ diff + 
             distanceKM_log +  
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m5 <- felm(count_emig2_log ~ diff + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m6 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)


### Output
stargazer(m1,m2,m3,m4,m5,m6,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          column.labels = c("Federal turnout"),
          column.separate = 6,
          dep.var.labels.include = T,
          dep.var.labels  = c("Emigration (log)"),
          covariate.labels=c("Difference",
                             "Distance (log kms)",
                             "Diff. manufacturing",
                             "Diff. services",
                             "Diff. GDP p.c.",
                             "Diff. UER",
                             "Diff. patents",
                             "Share manufacturing at origin",
                             "Share manufacturing at destination",
                             "Share services at origin",
                             "Share services at destination",
                             "GDP pc at origin",
                             "GDP pc at destination",
                             "UER at origin",
                             "UER at destination",    
                             "Patents pc at origin",
                             "Patents pc at destination",   
                             "GDP (log) at origin",
                             "GDP (log) at destination",
                             "Population (log) at origin",
                             "Population (log) at destination"),           
          digits = 3,
          add.lines = list(
            c("Origin-FE?", rep("YES", 6)),
            c("Destination-FE?", rep("YES", 6)),
            c("Year-FEs?", "YES", rep("YES", 6))),
          notes = c("Standard errors clustered by dyad"))




# Table A10: gravity model -- state election turnout ----------------------

### Regressions
grav$diff <- grav$diff_LT_turnout
m1 <- felm(count_emig2_log ~ diff + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m2 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m3 <- felm(count_emig2_log ~ diff + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m4 <- felm(count_emig2_log ~ diff + 
             distanceKM_log +  
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m5 <- felm(count_emig2_log ~ diff + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m6 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)


### Output
stargazer(m1,m2,m3,m4,m5,m6,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          column.labels = c("Federal turnout"),
          column.separate = 6,
          dep.var.labels.include = T,
          dep.var.labels  = c("Emigration (log)"),
          covariate.labels=c("Difference",
                             "Distance (log kms)",
                             "Diff. manufacturing",
                             "Diff. services",
                             "Diff. GDP p.c.",
                             "Diff. UER",
                             "Diff. patents",
                             "Share manufacturing at origin",
                             "Share manufacturing at destination",
                             "Share services at origin",
                             "Share services at destination",
                             "GDP pc at origin",
                             "GDP pc at destination",
                             "UER at origin",
                             "UER at destination",    
                             "Patents pc at origin",
                             "Patents pc at destination",   
                             "GDP (log) at origin",
                             "GDP (log) at destination",
                             "Population (log) at origin",
                             "Population (log) at destination"),           
          digits = 3,
          add.lines = list(
            c("Origin-FE?", rep("YES", 6)),
            c("Destination-FE?", rep("YES", 6)),
            c("Year-FEs?", "YES", rep("YES", 6))),
          notes = c("Standard errors clustered by dyad"))




# Table A11: gravity model -- state election result -----------------------

### Regressions
grav$diff <- grav$diff_LT_result
m1 <- felm(count_emig2_log ~ diff + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m2 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m3 <- felm(count_emig2_log ~ diff + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m4 <- felm(count_emig2_log ~ diff + 
             distanceKM_log +  
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m5 <- felm(count_emig2_log ~ diff + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m6 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)


### Output
stargazer(m1,m2,m3,m4,m5,m6,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          column.labels = c("Federal turnout"),
          column.separate = 6,
          dep.var.labels.include = T,
          dep.var.labels  = c("Emigration (log)"),
          covariate.labels=c("Difference",
                             "Distance (log kms)",
                             "Diff. manufacturing",
                             "Diff. services",
                             "Diff. GDP p.c.",
                             "Diff. UER",
                             "Diff. patents",
                             "Share manufacturing at origin",
                             "Share manufacturing at destination",
                             "Share services at origin",
                             "Share services at destination",
                             "GDP pc at origin",
                             "GDP pc at destination",
                             "UER at origin",
                             "UER at destination",    
                             "Patents pc at origin",
                             "Patents pc at destination",   
                             "GDP (log) at origin",
                             "GDP (log) at destination",
                             "Population (log) at origin",
                             "Population (log) at destination"),           
          digits = 3,
          add.lines = list(
            c("Origin-FE?", rep("YES", 6)),
            c("Destination-FE?", rep("YES", 6)),
            c("Year-FEs?", "YES", rep("YES", 6))),
          notes = c("Standard errors clustered by dyad"))




# Table A12: gravity model -- CSO density ---------------------------------

### Regressions
grav$diff <- grav$diff_CSO
m1 <- felm(count_emig2_log ~ diff + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m2 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m3 <- felm(count_emig2_log ~ diff + 
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m4 <- felm(count_emig2_log ~ diff + 
             distanceKM_log +  
             diff_manufacturing + diff_services + diff_gdppc + diff_UER + diff_patents +
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m5 <- felm(count_emig2_log ~ diff + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)
m6 <- felm(count_emig2_log ~ diff + 
             distanceKM_log + 
             orig_emp_VerarbGewerbe_share + dest_emp_VerarbGewerbe_share + 
             orig_emp_FinVersUnt_share + dest_emp_FinVersUnt_share + 
             I(orig_GDPpc/1000) + I(dest_GDPpc/1000) + 
             orig_unemprate_all + dest_unemprate_all + 
             orig_NPatents_erfinder_pc + dest_NPatents_erfinder_pc + 
             orig_GDP_log + dest_GDP_log + orig_pop_total_log + dest_pop_total_log | 
             orig_kreis + dest_kreis + year | 0 | dyad_FE2, data=grav)


### Output
stargazer(m1,m2,m3,m4,m5,m6,
          omit.stat = c("rsq", "f", "ser"),
          omit=c("AGS", "year"),
          star.cutoffs = c(0.05,0.01,0.001),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          column.labels = c("Federal turnout"),
          column.separate = 6,
          dep.var.labels.include = T,
          dep.var.labels  = c("Emigration (log)"),
          covariate.labels=c("Difference",
                             "Distance (log kms)",
                             "Diff. manufacturing",
                             "Diff. services",
                             "Diff. GDP p.c.",
                             "Diff. UER",
                             "Diff. patents",
                             "Share manufacturing at origin",
                             "Share manufacturing at destination",
                             "Share services at origin",
                             "Share services at destination",
                             "GDP pc at origin",
                             "GDP pc at destination",
                             "UER at origin",
                             "UER at destination",    
                             "Patents pc at origin",
                             "Patents pc at destination",   
                             "GDP (log) at origin",
                             "GDP (log) at destination",
                             "Population (log) at origin",
                             "Population (log) at destination"),           
          digits = 3,
          add.lines = list(
            c("Origin-FE?", rep("YES", 6)),
            c("Destination-FE?", rep("YES", 6)),
            c("Year-FEs?", "YES", rep("YES", 6))),
          notes = c("Standard errors clustered by dyad"))




# Figure A5: intraclass correlation coefficients --------------------------

### output container
out <- data.frame(matrix(nrow = 20, ncol=4))
colnames(out) <- c("sample", "outcome", "var.kreis", "var.resid")
out$sample <- factor(rep(1:4, each=5), levels=1:4,
                     labels = c("All municipalities", "Population < 1,000", "Population > 5,000", "Rural districts\n(Landkreise)"))
out$outcome <- factor(rep(1:5,4), levels=1:5,
                      labels = c("Federal election:\nturnout", "Federal election:\nLeft-right orientation", 
                                 "State election:\nturnout", "State election:\nLeft-right orientation",
                                 "CSO\ndensity"))


### regressions for all municipalities
# regressions
m.BT.turnout <- lmer(BT_turnout ~ 1 + (1 | AGSkreis), data = dmun)
m.BT.leftright <- lmer(BT_second_leftright_CHES ~ 1 + (1 | AGSkreis), data = dmun)
m.LT.turnout <- lmer(LT_turnout_lag ~ 1 + (1 | AGSkreis), data = dmun)
m.LT.leftright <- lmer(LT_votes_leftright_CHES ~ 1 + (1 | AGSkreis), data = dmun)
m.CSO <- lmer(CSO_stock_count_pc ~ 1 + (1 | AGSkreis), data = dmun)

# extract relevant variances
out[1,3] <- as.data.frame(VarCorr(m.BT.turnout))[1,4]
out[1,4] <- as.data.frame(VarCorr(m.BT.turnout))[2,4]

out[2,3] <- as.data.frame(VarCorr(m.BT.leftright))[1,4]
out[2,4] <- as.data.frame(VarCorr(m.BT.leftright))[2,4]

out[3,3] <- as.data.frame(VarCorr(m.LT.turnout))[1,4]
out[3,4] <- as.data.frame(VarCorr(m.LT.turnout))[2,4]

out[4,3] <- as.data.frame(VarCorr(m.LT.leftright))[1,4]
out[4,4] <- as.data.frame(VarCorr(m.LT.leftright))[2,4]

out[5,3] <- as.data.frame(VarCorr(m.CSO))[1,4]
out[5,4] <- as.data.frame(VarCorr(m.CSO))[2,4]



### regressions for municipalities below 5,000 inhabitants
# regressions
m.BT.turnout <- lmer(BT_turnout ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total < 1000)
m.BT.leftright <- lmer(BT_second_leftright_CHES ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total < 1000)
m.LT.turnout <- lmer(LT_turnout_lag ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total < 1000)
m.LT.leftright <- lmer(LT_votes_leftright_CHES ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total < 1000)
m.CSO <- lmer(CSO_stock_count_pc ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total < 1000)


# extract relevant variances
out[6,3] <- as.data.frame(VarCorr(m.BT.turnout))[1,4]
out[6,4] <- as.data.frame(VarCorr(m.BT.turnout))[2,4]

out[7,3] <- as.data.frame(VarCorr(m.BT.leftright))[1,4]
out[7,4] <- as.data.frame(VarCorr(m.BT.leftright))[2,4]

out[8,3] <- as.data.frame(VarCorr(m.LT.turnout))[1,4]
out[8,4] <- as.data.frame(VarCorr(m.LT.turnout))[2,4]

out[9,3] <- as.data.frame(VarCorr(m.LT.leftright))[1,4]
out[9,4] <- as.data.frame(VarCorr(m.LT.leftright))[2,4]

out[10,3] <- as.data.frame(VarCorr(m.CSO))[1,4]
out[10,4] <- as.data.frame(VarCorr(m.CSO))[2,4]



### regressions for municipalities above 5,000 inhabitants
# regressions
m.BT.turnout <- lmer(BT_turnout ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total > 5000)
m.BT.leftright <- lmer(BT_second_leftright_CHES ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total > 5000)
m.LT.turnout <- lmer(LT_turnout_lag ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total > 5000)
m.LT.leftright <- lmer(LT_votes_leftright_CHES ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total > 5000)
m.CSO <- lmer(CSO_stock_count_pc ~ 1 + (1 | AGSkreis), data = dmun, subset=pop_total > 5000)


# extract relevant variances
out[11,3] <- as.data.frame(VarCorr(m.BT.turnout))[1,4]
out[11,4] <- as.data.frame(VarCorr(m.BT.turnout))[2,4]

out[12,3] <- as.data.frame(VarCorr(m.BT.leftright))[1,4]
out[12,4] <- as.data.frame(VarCorr(m.BT.leftright))[2,4]

out[13,3] <- as.data.frame(VarCorr(m.LT.turnout))[1,4]
out[13,4] <- as.data.frame(VarCorr(m.LT.turnout))[2,4]

out[14,3] <- as.data.frame(VarCorr(m.LT.leftright))[1,4]
out[14,4] <- as.data.frame(VarCorr(m.LT.leftright))[2,4]

out[15,3] <- as.data.frame(VarCorr(m.CSO))[1,4]
out[15,4] <- as.data.frame(VarCorr(m.CSO))[2,4]


### regressions dropping Stadtkreise
# regressions
m.BT.turnout <- lmer(BT_turnout ~ 1 + (1 | AGSkreis), data = dmun, subset=nmun > 1)
m.BT.leftright <- lmer(BT_second_leftright_CHES ~ 1 + (1 | AGSkreis), data = dmun, subset=nmun > 1)
m.LT.turnout <- lmer(LT_turnout_lag ~ 1 + (1 | AGSkreis), data = dmun, subset=nmun > 1)
m.LT.leftright <- lmer(LT_votes_leftright_CHES ~ 1 + (1 | AGSkreis), data = dmun, subset=nmun > 1)
m.CSO <- lmer(CSO_stock_count_pc ~ 1 + (1 | AGSkreis), data = dmun, subset=nmun > 1)

# extract relevant variances
out[16,3] <- as.data.frame(VarCorr(m.BT.turnout))[1,4]
out[16,4] <- as.data.frame(VarCorr(m.BT.turnout))[2,4]

out[17,3] <- as.data.frame(VarCorr(m.BT.leftright))[1,4]
out[17,4] <- as.data.frame(VarCorr(m.BT.leftright))[2,4]

out[18,3] <- as.data.frame(VarCorr(m.LT.turnout))[1,4]
out[18,4] <- as.data.frame(VarCorr(m.LT.turnout))[2,4]

out[19,3] <- as.data.frame(VarCorr(m.LT.leftright))[1,4]
out[19,4] <- as.data.frame(VarCorr(m.LT.leftright))[2,4]

out[20,3] <- as.data.frame(VarCorr(m.CSO))[1,4]
out[20,4] <- as.data.frame(VarCorr(m.CSO))[2,4]


### compute ICC
out$var.total <- out$var.kreis + out$var.resid
out$ICC <- out$var.kreis / out$var.total


### Plot
p <- ggplot() + 
  geom_bar(data=out, aes(x=outcome, y=ICC*100, fill=outcome, group=1),
           stat="identity") + 
  geom_hline(yintercept = 50, lty=2) +  
  scale_x_discrete(limits=rev) +
  scale_y_continuous(limits=c(0,100), breaks=seq(0,100,20)) +
  scale_fill_manual(values = c("gray30","gray70", "burlywood3", "burlywood1", "firebrick3"), name="") +
  xlab("") +
  ylab("Intraclass Correlation Coefficient (ICC)") +
  facet_wrap(~sample, ncol=4) +
  coord_flip() +
  theme_bw() + 
  theme(strip.text = element_text(size=14, face="bold"),
        axis.title = element_text(size=14, face="bold"),
        #axis.text.x = element_text(size=12, angle=90, hjust=1, vjust=0.5),
        axis.text = element_text(size=12),
        legend.position = "none") 
p
ggsave("FigureA5_ICC_mun_county.png", p, width=12, height=5)





# Table A13: Little evidence of contextual effects ------------------------

### See STATA do-file


